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Abstract. We study nonlinear waves in Newton's cradle, a classical mechanical 
system consisting of a chain of beads attached to linear pendula and interacting 
nonlinearly via Hertz's contact forces. We formally derive a spatially discrete 
modulation equation, for small amplitude nonlinear waves consisting of slow mod- 
ulations of time-periodic linear oscillations. The fully-nonlinear and unilateral 
interactions between beads yield a nonstandard modulation equation that we call 
the discrete p-Schrodinger (DpS) equation. It consists of a spatial discretization 
of a generalized Schrodinger equation with p-Laplacian, with fractional p > 2 de- 
pending on the exponent of Hertz's contact force. We show that the DpS equation 
admits explicit periodic travelling wave solutions, and numerically find a plethora 
of standing wave solutions given by the orbits of a discrete map, in particular spa- 
tially localized breather solutions. Using a modified Lyapunov-Schmidt technique, 
we prove the existence of exact periodic travelling waves in the chain of beads, 
close to the small amplitude modulated waves given by the DpS equation. Us- 
ing numerical simulations, we show that the DpS equation captures several other 
important features of the dynamics in the weakly nonlinear regime, namely modu- 
lational instabilities, the existence of static and travelling breathers, and repulsive 
or attractive interactions of these localized structures. 

1. Introduction and main results 

This paper concerns the mathematical analysis and numerical simulation of non- 
linear waves in the Newton's cradle, a classical mechanical system consisting of a 
chain of beads suspended from a bar by inelastic strings (see figured]). All beads 
are identical and behave like a linear pendula in the absence of contact with nearest 
neighbours, but mechanical constraints between touching beads introduce geometric 
nonlinearities. 

In the last decades, this system has been considered as a reference problem to 
test multiple impacts laws, aimed at evaluating e.g. the post-impact velocities of 
all beads after a bead is released at one end of the cradle [131 121 E2J [53] (see also 
[40J for additional references). This problem is very delicate, because the collisional 
dynamics involves nonlinear elastic waves that propagate along the granular chain. 
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One of the important factors that influence wave propagation and multiple colli- 
sions is the nature of elastic interactions between beads (see [361 Ell EH] and refer- 
ences therein). Hertz's theory [5HH6] allows to compute the repulsive force between 
two identical and initially tangent spherical beads that are compressed and slightly 
flatten. When the distance between their centers decreases by S ~ (see figured]), 
the repulsive force / is f(S) = k5 a at leading order in 5, where k depends on the 
ball radius and material properties and a = 3/2. This result remains valid for much 
more general geometries [T2"| HB] . but a can be larger for irregular contacts (a = 2 
in the presence of conical asperities [31]) or smaller for surfaces that squeeze more 
easily (a = 1 for solid cylinders in contact at two ends). Hertz's type contact forces 
have several properties that make the analysis of wave propagation particularly dif- 
ficult. Firstly, they consist of unilateral constraints i.e. no force is present when 
beads are not in contact. Moreover, they are fully nonlinear for a > 1, and in that 
case classical linear wave theory becomes useless. In addition f"(0) is not defined 
for a < 2, therefore the use of perturbative methods is more delicate. Developping 
analytical tools to overcome these obstacles is important, because the latter are at 
the heart of wave propagation in granular media. 

5 

Figure 1. Left : Newton's cradle. Right : schematic representation 
of two compressed beads. 



A simplified model for Newton's cradle reads in dimensionless form 

cPx 

+ x n = V'(x n+1 - x n ) - V'(x n - x n -i), neZ, (1) 

where x n (t) G 1R is the horizontal displacement of the nth bead from the ground- 
state, i.e. the equilibrium position at which each pendulum is vertical, beads are in 
contact at a single point and uncompressed. The interaction potential V takes the 
following form as x — > 

V(x) = V (x)+H(-x)W(x), V (x) = —— \x\ 1+a H(-x), W (x) = o(\x\ 1+a ) , 

1 + a 

(2) 

where H denotes the Heaviside function vanishing on R and equal to unity on 
M + , o is the usual Landau's symbol and a > 1 a fixed constant. The potential Vq 
corresponds to a generalized Hertz contact force and W incorporates higher-order 
corrections. System ([2]) is Hamiltonian with total energy 

n = J2 1 2^y '' + l< + V(x n+1 -x n ). (3) 
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Equation ([T]) considered on the infinite lattice Z applies in principle to large ensem- 
bles of beads (so that boundary effects can be neglected for waves propagating in 
their core), but this model is also suitable to describe strongly localized standing 
waves in relatively small systems. 

The simplest model to encompass difficulties inherent to Hertz's contact nonlin- 
earities consists of a chain of identical spherical beads without local oscillators, in 
contact with their neighbours at a single point when the chain is at rest. Intensive 
research has been carried out on the analytical description of compression pulses 
in this system. Nesterenko analyzed the problem using a formal continuum limit 
[58| [59] . and found an approximate pulse solution consisting of a compression solitary 
wave with compact support (see also [3] for a slightly different continuum limit, and 
[63J for an extension to dimer chains). In addition, Chatterjee [2] has subsequently 
obtained an improved approximation revealing a double-exponential decay of the 
wave profile. As shown by MacKay [56] (see also [15]), exact solitary waves exist 
since a theorem of Friesecke and Wattis [30] proving their existence readily applies 
to the chain of beads with Hertz contact forces. In addition, English and Pego [19] 
have shown that these solitary waves have a doubly-exponential decay in the case of 
Hertz's force. More details on the dynamics of these solitary waves can be found in 
the review [62] • In addition, much more properties of solitary waves are known when 
an external load /o is applied at both ends of the chain and all beads undergo a small 
compression 8 . The dynamics around this new equilibrium state reduces to the one 
of the classical Fermi- Pasta-Ulam (FPU) lattice [59, 69J, in the case of which soli- 
tary waves with exponential decay [301 ESI [26], EH STJ 021 EH E7J 135] and two-soliton 
solutions [37] are known to exist, small amplitude solitary waves are well described 
by the KdV equation [261 El S3 EH UJ and are stable [271 M ESS EH EH [39] . 

This mapping between the chain of compressed beads and the FPU model is 
interesting, because the latter is known to display a rich dynamical behaviour and 
sustain many other kinds of nonlinear waves (see [HI [32] for recent reviews on this 
topic). Among the most fundamental excitations existing in FPU chains, periodic 
travelling waves [201 HD Ell EH [35] are particularly important to understand energy 
propagation and dispersive shocks [TTJ. However, no solutions of this type are known 
for the the chain of beads in the absence of external load. The above mentioned 
periodic travelling waves are degenerate when f — > 0, because the sound velocity 

1 /6 

(i.e. the maximal velocity of linear waves) vanishes as f . For this reason the 
uncompressed chain of beads is commonly denoted as a "sonic vacuum" [59] . 

In contrast with the above statement, we prove in this paper that nonlinear pe- 
riodic travelling waves exist in Newton's cradle, due to the interplay between the 
on-site oscillators and fully-nonlinear contact interactions among beads. 

Theorem 1. Consider a potential V G C 2 (M.) taking the form (T^) with a > 1 and 
W"(x) = o(|x| a_1 ) as x 0. There exists a > such that for all a G (0, a ) and 
q G (— 7r, 7r], system admits a periodic travelling wave solution 




(4) 
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with amplitude a and wavenumber q, where the wave frequency u a > 1 satisfies a 
nonlinear dispersion relation 

w. = 1 + - a*' 1 I sin (|)r +1 + o(a"- 1 ), (5) 

v^( a 2 -i)r(^) 

T = _n„TV^ > ( 6 ) 



o;2 a r(f) 



is 



and r(x) = J °° e 1 d£ denotes Euler's Gamma function. The function v a 
2n -periodic, odd and belongs to C 2 (M). It takes the form 

v a (0 = K h P h [V '(sm (£ + q) - sin£) - V£(sin£ - sin (f - g))] + R a ($, (7) 

where \\Ra\\L°°(R) ~ as a — > and the linear operators Ph, Kh are defined by 

(P h f)(0 = f(0-(^0- r sin (s)f(s)ds, (8) 

7T ./n 
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(AT h /)(0 = (sinO - / (s - tt) cos (s) /(s) (is + / sin (f - s) /(s) ds. (9) 

These waves display several unusual features. Firstly, the family of periodic trav- 
elling waves (parametrized by q and a w 0) is singular when a 6 (1,2) (which is 
the case for Hertz's contact law), in the sense that the frequency u a defined by (jSJ) 
is not differentiable with respect to a at a = 0. This originates from the limited 
smoothness of the potential Vq at the origin [V (CT) is not defined). Moreover, 
as a result of fully nonlinear interactions between beads, for all wavenumber q the 
wave frequency u a converges towards unity in the small amplitude limit a — > (i.e. 
the linear phonon band reduces to a single frequency). 

The periodic travelling waves (H])-© carry an energy flux in the direction of 
wave propagation when q £ {0,tc}. In the special case q = 0, solution (j4]) re- 
duces to x n (t) = —a sin (t) and corresponds to in-phase oscillations of uninteract- 
ing linear pendula. For q = tt these solutions correspond to binary oscillations 
x n (t) = a (-1)" +1 (sin (cj a t) + a a - l v a {u a t)). 

The proof of theorem [1] proceeds in two steps. Firstly, we formally derive in 
section 12.11 an amplitude equation 

3A 

2^To— - = (A n+1 - An) \A n+1 - A n \ a ~ l - (A n - A n ^) \A n - A n ^\ a ~\ (10) 
or 

which describes small amplitude approximate solutions of ([T]) taking the form 

x n pp (t) = e(A n (r)e lt + A n (r)e- u ), r = e a ~\ (11) 

where e > is a small parameter and A n (r) G C. The approximate solutions 
determined by ( !IU|) - (fTT|) consist of slow modulations in time of solutions of problem 
(JTJ linearized at x n = 0. The ansatz ( ITT]) and amplitude equation (TTTJj) are consistent 
with the nonlinear problem (jTJ, in the sense that ( II ip approximately satisfies (|T]) up 
to an o(e a ) error term. We show in section 12.21 that (llOp admits an explicit family 
of periodic travelling wave solutions, which provide the harmonic part of (jlj) and 
the leading order terms of the dispersion relation fl5]). In a second step, we prove in 
section |3]the existence of exact periodic travelling wave solutions of (JTJ close to these 
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approximate solutions. For this purpose we consider an advance- delay differential 
equation that determines travelling waves of given period. We solve this problem by 
adapting the Lyapunov-Schmidt technique to the kind of nonlinearity with limited 
smoothness present in ([!]) (using the contraction mapping theorem instead of the 
implicit function theorem employed in the usual case). 

Our approach has the advantage of giving explicitly the principal part of travelling 
waves in the small amplitude limit, since expressions (JI])-(J7]) provide an approxima- 
tion of the exact solutions up to an o(a a ) error term. In fact, we will show that far 
more complex approximations of order o(a^) can be obtained for all (3 > 0, which 
constitutes an unexpected result given the limited smoothness of V at the origin. 

In section 14.11 we complete theorem [I] by numerical results, in order to study the 
limits of the above analysis when one leaves the small amplitude regime. Along the 
same line, it would be interesting to analyze periodic travelling wave solutions of 
([T]) using variational methods j2Ul EES! ED ES] or degree theory [20] , since this should 
relax the assumption of small amplitude waves (with the limitation of describing 
travelling wave profiles and dispersion relations with less precision than perturbative 
methods). 

We call equation (JTQ|) the time- dependent discrete p-Schro 'dinger (DpS) equation, 
because it can be seen as a finite-difference spatial discretization of the generalized 
Schrodinger equation i d T A = A p A, where the usual Laplacian operator is replaced 
by the one- dimensional p-Laplacian A P A = d^A \d^A\ p ~ 2 ) and p = a + 1. Its 
fully nonlinear structure and fractional power nonlinearities originate from the fully 
nonlinear interaction forces of ([1]). Interestingly, the unilateral character of Hertz's 
contact forces is averaged out along the fast oscillations of the pendula, which results 
in much simpler nonlinearities at the level of the DpS equation. In the present 
context, the DpS equation substitutes to the classical discrete nonlinear Schrodinger 
(DNLS) equation 

dA 

«^ = A n+1 -2A n + A n _ 1 ±A n \A n \ 2 , (12) 
or 

which can be derived for Hamiltonian lattices with smooth (at least C 4 ) potential 
energies and linear dispersion (HI [151 EO HB] • 

Equation ffTUl) tells in fact much more on the dynamics of Newton's cradle than 
what is described in theorem [TJ In section 14.21 we numerically observe that certain 
periodic travelling waves of ([[]) are unstable through the phenomenon of modula- 
tional instability. Their envelope self-localizes under the effect of small perturba- 
tions, yielding spatially localized and time-periodic intermittent compressions of the 
beads. These localized oscillations propagate along the chain and sometimes remain 
pinned at some lattice sites. These waves correspond to discrete breathers, a class of 
nonlinear excitations ubiquitous in spatially discrete systems [HJ E5J H21 El HU EE]- 
They have been studied both experimentally and theoretically in different types of 
granular chains, in the absence of local oscillators and under precompression [XU1I73] . 
We check numerically that the DpS equation accurately describes the modulational 
instability of certain small amplitude periodic travelling waves, provided one avoids 
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near-critical cases and higher-harmonic instabilities (section |4.2p . Moreover, numer- 
ical computations reveal that breather solutions of the DpS equation (determined 
as homoclinic orbits of a discrete map) describe with high precision the envelope 
of small amplitude breathers of ([1]) (sections 12.21 and 14.31) . The DpS equation also 
qualitatively reproduces differents kinds of slow interactions of discrete breathers, 
namely the fissionning of a localized perturbation into slowly travelling breathers, 
and the merging of two breathers into a single one after a long transcient (section 
14. 3p . More generally, we numerically find a plethora of standing wave solutions of the 
DpS equation given by the orbits of an area-preserving and reversible map (sections 
12. 2p . which suggests the existence of many time-periodic standing wave solutions of 
(PQ) with a huge variety of spatial behaviours. 

The paper is organized as follows. Section [2] contains the formal derivation of 
the DpS equation and a short discussion of some travelling wave and standing wave 
solutions. The proof of theorem [1] is given in section [3] and most numerical com- 
putations performed in section HJ Lastly, section [5] discusses new perspectives and 
open problems resulting from the present work. 

2. The DpS equation for modulated waves 

The aim of this section is twofold. Section 12.11 presents a formal derivation of 
the DpS equation, starting from Newton's cradle equations ([I]). One has to stress 
that this gives by no means a justification that the DpS equation approximates the 
original system for appropriate initial conditions and timescales (these problems will 
be studied analytically in section 13.11 for exact travelling waves, and numerically in 
section H] for various initial conditions). In section fl?2\ we study particular classes 
of solutions of the DpS equation, i.e. periodic travelling waves and standing waves. 
Periodic travelling waves are computed explicitly, a result which will serve as a 
basis for the analytical and numerical studies of sections [3] and HI Standing waves 
are studied numerically for a = 3/2, with a special emphazis on spatially localized 
standing waves (these solutions are important to describe modulational instabilities 
in system ([1]) , as section H] will show) . 

2.1. Formal derivation and elementary properties. In this section we formally 
derive the DpS equation (flQj) from the original dynamical equations ([1]), using a 
multiscale expansion technique. We look for solutions of ([T]) in the form of slow 
modulations in time of 27r-periodic functions (2tt being the period of the linear 
on-site oscillators). More precisely, we set 

x n {t)=X n (r,t), (13) 

where r = fit corresponds to a slow time, fi > is a small parameter, and 

X n (r,t + 2n) =X n (r,t) (14) 

for all t, t 6 M. We replace flTJ by the infinite system of coupled PDE 

[(jid T + d t ) 2 + l}X n = V f (X n+1 -X n )-V , (X n -X n _ 1 ), n6Z, (r,i)6l 2 (15) 

with periodic boundary conditions (j!4p . The restriction to the line r = fit of any 
solution of ffl5l) defines a solution of ([1]). Now we look for a family of small amplitude 
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solutions of ffl5|) corresponding to nearly harmonic oscillations 

X%t, t) = e (Al(r) e u + A n {r) e~ u ) + e^(r, t) (16) 

and parametrized by e > close to 0. The complex amplitude A e n varies slowly in 
time, but no assumptions are made on its spatial behaviour. The remainder R n e K 
satisfies f * R n (r,t) e ±lt dt = and we fix (3 > 1. The assumptions made on V in 
theorem [I] imply 

V\x) = V Q \x) + H{-x)o{\x\ a ), V \x) = -\x\ a H(-x). (17) 

This leads us to fix fl = e"" 1 and j3 = a, so that in equation ( fl5l) the time-modulation 
term jld 2 t X n , the nonlinear interaction forces V and the remainder e^R e n have the 
same order e a . Setting A e n = A n + o(l) and R e n = R n + o(l) as e — > 0, inserting ( FIBj) 
in equation ( fl5l) and multiplying the latter by e~ a , one obtains 

2id T A n e u - 2id T A n e~ lt + (d 2 + l)R n (18) 
= V '[ (A n+1 - A n ) e u + c.c. } - Y Q [ (A n - A n ^) e u + c.c. ] + o(l) 

as e — )■ 0, where we denote by c.c. the complex conjugates. Now assume that a 
family of solutions X e n exists for all e ~ 0. Letting e — > in the above equation 
yields 

2icU n = /(A n+1 - A n ) - f(A n - A n _0, (19) 

f(z) = — / Vft( ze' + z e" 4 * ) e" rf 
2tt J 

by projection on the first Fourier harmonic. Moreover, for all fixed r the 27r-periodic 
function R n {r, .) is determined as a function of A n ±i{j) and A n (r), and defined as 
the unique 27r-periodic solution of 

( d 2 t + 1 ) R n = V{[ - A n ) e if + c.c. ] - V£[ {An - A n _{) e u + c.c. ] - 2id T A n e u + c.c. 

being L 2 -orthogonal to e ±lt . Now there remains to evaluate f(z). Setting z = re td 
and using the change of variable s = t + 9 in the integral defining /, one obtains 



f(*) = i- I V ! (2rcoss)e- is ds 
2tt J t 



where T denotes any interval of length 2n. Now we fix T = (— tt, 7r) and use the fact 
that Vq(x) = —\x\ a H(—x). One obtains after elementary computations 

f(z) = 2 a - 1 c a z\z\ a - 1 (20) 



where 



o 



c a = - [ 1 (cost) a+1 dt (21) 



(22) 



7T 

is a Wallis integral with fractional index a + 1. It follows that 

ir(|)r(f + i) 



Q+l 
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where T(x) = / +o ° e t t x 1 dt denotes Euler's Gamma function (see [T\, formula 6.2.1 
and 6.2.2 p. 258). Since T(l/2) = y/n and T(a + 1) = aT(a), we obtain finally 



2aT{ 



a ' 



^ A(^-i)r(^)' (23) 

As a conclusion, introducing r = 2 1 ~ Q c~ 1 , equations (|T9|) -( l20|) yield 
dA 

2iT Q —^ = {A n+l -A n )\A n+1 -A n \ a - 1 -{A n -A n _ 1 )\A n -A n _ 1 \ a -\ n e Z, (24) 

OT 

i.e. one recovers the DpS equation ( TTUj) introduced in section [1] and expression ([6]) 
of coefficient r . 

Note that in the particular case of Hertz's law, one has 

a = 3/2, r = « 1.545 

(reference pQ, formula 6.1.32, 6.1.10 p. 255-256). Moreover, in the case of nearly 
linear unilateral interaction forces (i.e. when a — > 1 + ), one can use expression ( 12T|) 
to obtain lim tq = 2. 



Remark 1. The discrete quasilinear Schrddinger (D-QLS) equation 
.dA n 

was derived in reference [64] for Hamiltonian lattices with 



k [ {A n +i ~ A n ) \ A n+ i — A n \ 2 — (A n — A n _i) \ A n — A n _i\ 2 ] + A n |v4„| 2 (25) 



U = £ \ ^ + \ X l-\< + k (*n + l - *uY (26) 

(k > 0), where the purely quartic interaction potential yields a nonlinear dispersion 
and the on-site potential is anharmonic. For stationary solutions A n (r) = a n e inr , 
the additional local nonlinearity makes equation l[25\) easier to analyze than (24\ ), 
because can be analytically studied near the anticontinuum limit k — > [51 H]. 
In addition, the existence of discrete breather solutions of ( f^5]j has been proved in 
reference [64], as well as several results concerning their stability, spatial decay and 
continuum limit. Unfortunately, system < f^3j) does not apply to the Newton's craddle 
problem, because nonlinear terms of the local potential U(x) = 1 — cos x of nonlinear 
pendula are of higher order than the classical Hertz contact interactions (for this 
reason on-site nonlinearities have been neglected in Hamiltonian $2$). 

Equation (124)1 possesses some interesting properties, in particular conserved quan- 
tities. Let us consider spatially localized solutions of ( 124|) satisfying {A n (r)} e ^(Z). 
One can check that the squared I2 norm Ylnez l^«| 2 ^ s invariant by the flow of fl24|) . 
In addition, ( I24p can be written as a Hamiltonian system 

a 4 _ 9H_ _dH_ 

OPn oA n 
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with p n = i A n and the Hamiltonian 

h = -- — — t— lAi+i - A n \ a+i . 

(a + l)r ^ 

It follows that the £ a +i norm of the forward difference {A n+ \— A n } is also a conserved 
quantity. If {A n (r)} £ £i(Z) then there exists a third conserved quantity 



In addition the DpS equation admits the gauge invariance A n — > A n e llfi , the trans- 
lational invariance A n — >■ v4 n + c and a scale invariance, since any solution A n (r) 
of f )24|) generates a one-parameter family of solutions a v4 n (|a| a_1 r), a £ JR.. This 
property explains the structure of velocity-amplitude relations for travelling wave 
solutions and the frequency-amplitude scaling of standing wave solutions that will 
be obtained in the next section. 

To end this section, we point out an interesting formal continuum limit of (|2"3|) 
for solutions A n (r) varying slowly in space and time with an appropriate scaling. 
Setting p = a+ 1, A n (r) = s) with f = hn, s = (2r Q )- 1 h a+l T in equation flU} 
and letting h go to 0, one obtains the generalized Schrodinger equation 

id s ip = A p ?p, £ £ R, (27) 

where the usual Laplacian operator is replaced by the one-dimensional p-Laplacian 
A p ip = d%( d^ip \d^ip\ p ~ 2 ). In addition, more general fully nonlinear generalized 
Schrodinger equations would result from different continuum limits, setting A n (r) = 
B(£,t) e iqn . 

2.2. Time-periodic travelling and standing waves. It is well known (see e.g. 
|18j ) that the classical DNLS equation f fT2|) admits explicit periodic travelling wave 
solutions of the form 

A n (r) = R e i(nr- q n-<t>)^ R > ^ ^ 

where the frequency Q is fixed by the amplitude R and the wavenumber q through 
a nonlinear dispersion relation. The same property holds true for the DpS equation 
QUI) due to its gauge invariance. Inserting the ansatz ([2"5]) in (fTOl) . one obtains after 
elementary computations 

T tt = 2 a i? Q - 1 |sin(|)| a+1 . (29) 

Since (TTDT) was derived as an amplitude equation from the original system flTJ, com- 
bining the ansatz (j28p with f ll3p -f ll6p yields approximate periodic travelling wave 
solutions of (Cp). These solutions take the form 

x%(t) = e(A n (e a -H)e u + A n (e a -H)e' u ) (30) 

= a cos (qn — u tw t + 0), (31) 

where we denote by a = 2eR > the wave amplitude and by 

u tw = l + -a a - 1 \sm(l)\ a+1 (32) 

its frequency. 
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Interestingly, these expressions are exact for q = and correspond to linear in- 
phase oscillations of all oscillators of ([I]) at frequency u tw = 1. More generally, as 
announced in theorem [1] there exist in fact exact periodic travelling wave solutions 
of ([1]) close to the above approximate solutions. The proof of this statement will be 
the object of section [3j 

The nonlinear dispersion relation displays an unusual feature, since for all 
wavenumber q the wave frequency u tw converges towards 1 in the small amplitude 
limit a —7- 0. This is due to the fact that dispersion originates from fully nonlinear 
interactions between beads. Additionally, one can notice that the band of frequencies 
u tw of these periodic travelling waves lies above u = 1. 

Another special class of solutions of ( TTUj) takes the form 

A n (T)=R n e i( fi T+ *\ R n ER. (33) 

Introducing a n = 2eR n and u sw = 1 + fie a_1 , the approximate solutions deduced 
from ( Tl6|) read 

4 w (t) = a n cos(w sw t + 0), (34) 
and the DpS equation (fTOj) becomes 

— ix a n = (a n+ i — a n ) \a n +i — a n \ a 1 — (a n — a n _i) \a n — a n _i| Q , (35) 

with 

H = (uj sw - 1) 2 a r . 

Solutions of the form (|34p correspond to standing waves, i.e. they oscillate periodi- 
cally with time and their nodes and extrema do not change. We shall refer to ( 135)) 
as the (real) stationary DpS equation. 

Equation fl35|) has been previously derived for special classes of lattices that sus- 
tain standing wave solutions x n (t) = a n <p(t) with space-time separation (see re- 
mark |5] for more details). 

Remark 2. Equation / T53]) can be simply derived in some class of nonlinear lat- 
tices with special interaction potentials. Consider the symmetrized potential U(x) = 
Vq{— \x\) = (1 + a) -1 |s| a+1 satisfying U'(x) = x and the nonlinear system 

d?x 

—^ + Kx n = U'(x n+1 -x n )-U'(x n -x n -.i), neZ. (36) 
at z 

System IfSh)) admits some simple mechanical interpretations. When k > 0, it de- 
scribes e.g. the oscillations a chain of linear pendula coupled by anharmonic torsion 
springs with potential U . For k < the same analogy can be carried out, replacing 
the linear pendula by inverted pendula. For k = 0, system IfSh)) reduces to the Fermi- 
Pasta-Ulam lattice with potential U, that represents e.g. a chain of masses coupled 
by anharmonic springs. Equilibria of IfSh}) satisfy equation If35*\) with fi = —k. More- 
over, due to the special form of U, one can find standing wave solutions of l[3b}) 
taking the form x n (t) = a n ip(t), where if is a periodic solution of 

-^+Kip + yU'{tp) =0 
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and {a n } is a bounded solution of j[35*\) . For k = and when a > 3 is an odd integer, 
spatially localized standing wave solutions of IfSh)) have been studied in the physics 
literature j50j [22j (2TJ [70] using equation l[35\) . The case k > and a = 3 /ias 6een 
studied in [23] a/ong t/te same lines. 

The simplest case of equation (135]) corresponds to // = (i.e. w sw = 1), where it 
follows that a n+ \ — a n is constant, and thus all solutions of (135]) read a n = An + (3, 
A, (3 G R. In that case, (134"]) yields in fact exact solutions of (JT|) taking the form 
x n {t) — (An + /3) cos (t + 0) and corresponding to collective in-phase oscillations. 

For n 7^ 0, it is interesting to note that a n = l/i] 1 ^ a n satisfies the renormalized 
equation 

- sign(/i) a n = (a n+ i - a n ) \a n+1 - an]"' 1 - (a n - a n _i) \a n - a n _i| Q_1 , (37) 

where by definition sign(/i) = 1 for fi > and sign(/i) = —1 for /i < 0. Consequently 
it is sufficient to restrict to the cases /i G { — 1,0,1} to find all solutions of fl3"5]) . In 
what follows we fix = ±1 and drop the tilde in equation ( )37j) for notational 
simplicity. 

Multiplying flHTj) by a n and summing over n, one obtains for all {a n } G ^(^) 

where an index change has been performed to simplify the right-hand side. Con- 
sequently, nontrivial localized solutions {a n } G can be found only for /i > 0, 
and the same conclusion holds true for periodic solutions (just slightly modifying 
the above computation). As a consequence we shall restrict ourselves to the case 
(j, = 1. 

In order to reformulate (1371) as a two-dimensional mapping, we introduce the 
auxiliary variable 

b n — { a n ~ a n-l) \ a n ~ a n-l| Q ■ (38) 

Inverting (138]) and shifting the index n, one finds 

a n+1 = a n + b n+ i |6„+i|~ _1 . (39) 
Moreover, equation ( )37j) at rank n + 1 reads 

b n +2 = b n+ i - a n+1 . (40) 

Introducing the variable 

U n = (a n , b n+ i), 

equations (I39l - (f40l) define a two-dimensional mapping U n+ \ = G{U n ). One can 
check that detDG{U) = 1 for all U G M 2 \ {0}, which implies that G is an area- 
preserving map. Moreover, the mapping possesses the invariance U n — > —U n since 
G(—U) = —G(U). In addition the map G is reversible with respect to the symmetry 
R defined by R(a, b) = (a, —a — b), i.e. one has RG~ Y = G o R. Equivalently, {U n } 
is an orbit of G if and only if {RU_ n } is an orbit of G. This property originates 
from the invariance n — > —n of ( )37j) . 
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In what follows we consider f l39]) - (pi0]) for a = 3/2. Several orbits of G are shown 
in figures [2] and [3j Figure [2] shows an orbit of G, for an initial condition with a® = 
and h\ very small. The trajectory suggests the existence of pairs of symmetric 
orbits homoclinic to 0, i.e. satisfying lim U n = 0. Figure [3] shows many kinds of 

n— »±oo 

trajectories frequently encountered in area-preserving or reversible mappings. Plots 
in the left column show periodic orbits and invariant tori, while right plots reveal 
intricate trajectories in the vicinity of stable and unstable manifolds. 




Figure 2. Orbit of G for the initial condition a = 0, b\ = 10 



Now we concentrate on solutions homoclinic to 0. To approximate these solutions 
numerically, we consider equation ( 1371) with periodic boundary conditions a n+ ^ = a n 
(computations are performed for iV = 50), and look for solutions close to homoclinic 
ones. This approach is consistent with figure 131 which shows periodic orbits and in- 
variant tori very close to the pair of homoclinics. We numerically solve this nonlinear 
equation using the fsolve function of the software package Scilab, based on a modified 
Powell hybrid method. 

We start by computing an homoclinic solution of ( I57)) with site-centered symmetry 
a N/2-n = O'N/2+n- The numerical iteration is initialized by setting a N / 2 = —a and 
a n = elsewhere. We fix a = 0.1, of the order of the size of the homoclinic loop of 
figure [3J In that case the iteration converges towards a spatially localized symmetric 
solution {a n }, whose profile is shown in the top left plot of figure HI As shown by 
the top right plot in semi- logarithmic scale, the homoclinic solution has a super- 
exponential spatial decay. This unusual feature comes from the fact that G is not 
differentiable at the fixed point U = (since a > 1), hence the classical results on 
exponential convergence along stable and unstable manifolds do not apply. 

One can also compute an homoclinic solution of (137]) with bond-centered symme- 
try a N / 2 -n = — a N/2+n+i- in that case one starts the numerical iteration by setting 
ajv/2 = a, cln/2+i = ~ a and a n = elsewhere. The iteration converges towards a 
spatially localized antisymmetric solution {a n } shown in figure H] (bottom left plot), 
which also decays super-exponentially. More complex homoclinic solutions can be 
computed by changing the initial condition (one example is shown in the bottom 
right plot) and the lattice size. 
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Figure 3. Trajectories of (|3l?]) - (l4"0"]) for a = 3/2, in neighbourhoods 
of U = of different sizes (top : ||Z7 n ||t» < 0.1, middle : ||?7n||oo < 1; 
bottom left : ||?7 n ||oo < 8, bottom right : ||f7 n ||oo < 100). Left plots 
shows periodic orbits and invariant tori, and right plots mainly focus 
on trajectories in the vicinity of stable and unstable manifolds. Marks 
in the first and second plot of the right column correspond to simple 
periodic orbits that can be explicitly computed, a period-2 orbit a n = 
2 T= « (— l) n , a period-3 orbit {a n } = {. . . , 0, a, —a, 0, a, —a, . . .} with 
a = (1 + 2 a ) T= ^ , a period-4 orbit {a n } = {. . . , 0, —a, 0, a, 0, —a, . . .} 
with a = 21^. The trajectory of figure [2] near the double homoclinics 
is visible in the top right plot, surrounding the period-2 orbit. 

Any solution of (I3TI) corresponds to the approximate solution of ([TJ 

x s ™(t) = an cos (u; sw t), fj, = (u sw - 1) 2%. (41) 

In particular, any orbit homoclinic to corresponds via (HTj) to a one-parameter 
family of time-periodic and spatially localized approximate solutions of (0Q), i.e. 
discrete breathers. Interestingly, it follows from expression (I4ip that the spatial 
extension of these breathers is independent of their amplitude (or equivalently of 
their frequency). This property is unusual (the spatial extension of classical discrete 
breathers diverges when their amplitude goes to |43J HI]) and originates from 
the fully-nonlinear interactions between beads present in ([1]). It is an analogue 
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Figure 4. This figure shows three different solutions of ( 1371) homo- 
clinic to 0. Top left plot : spatially symmetric solution. Top right 
plot : amplitude of the same homoclinic solution in semi- logarithmic 
scale, showing its super-exponential decay (profiles are only plotted for 
18 < n < 34, because \a n \ drops below machine precision at the other 
sites). Bottom left plot : spatially antisymmetric solution. Bottom 
right plot : double-humped homoclinic solution. 



for discrete breathers of a similar invariance present in Nesterenko's solitary wave 
[58"! I5T?] , whose spatial extension is independent of amplitude and velocity. 

More generally, this study leads us to conjecture the existence of time-periodic 
standing wave solutions of (QQ) having a large variety of spatial behaviours, e.g. 
spatially-periodic (with arbitrarily large periods), quasi-periodic, spatially localized, 
homoclinic to periodic or spatially disordered. 

Remark 3. Reference |40] provides interesting numerical and experimental results 
on standing waves in Newton's cradle with a small number of beads, in particular its 
relaxation towards in-phase oscillations when energy dissipation during collisions is 
not neglected, and the near-recurrence of specific modes of oscillation in the conser- 
vative case. 

3. Proof of theorem Q] 

The aim of this section is to prove theorem [Tj, i.e. show the existence of small 
amplitude exact periodic travelling wave solutions of ([T]) close to the approximate 
travelling waves provided by the DpS equation. For this purpose we reformulate the 
search of periodic travelling waves as a fixed point problem in section 13.11 using a 
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suitable scaling of the solutions deduced from the formal analysis of section [2j In 
section [372] we solve the resulting equation by the contraction mapping theorem. We 
end this section by pointing out some symmetry properties of periodic travelling 
waves (section [373]) . 

3.1. Fixed point problem for periodic travelling waves. We search for solu- 
tions of ([T]) taking the form x n (t) = w(£), where u is 27r-periodic, £ = qn — cut 
denotes a moving frame coordinate, and q G [— 7T, 7r], oj G R* are the wave parame- 
ters. Equation ([[]) yields the advance-delay differential equation 

u;V'(0 + u(0 = V'(u{£ + q)- u(0) - V'{u(C) - u(£ -q)), £ G R. (42) 

Notice that (142]) possesses the symmetry w(£) — > —u(—£). To simplify the analysis, 
we look for solutions of (142]) invariant under this symmetry, i.e. odd in £. 

We look for solutions with frequency u close to 1 as in section 12.21 Problem ( 142]) 
can be rewritten in more compact form 

u" + u=(l-w 2 )u" + N(u), (43) 

where 

[#(«)] (0 = V'Mt + Q) ~ - V'(u{0 - u{£ - q)). (44) 

Now we reformulate equation (143]) as a suitable fixed-point problem for nontrivial 
solutions (i.e. our formulation will eliminate the trivial solution u — 0). 
We begin by introducing some notations. We consider the Banach space 

X = { W GL7 p ° er (0,27r), «(-£) = } 

endowed with the L°° norm, where Cp er (0, 2ir) denotes the classical Banach space of 
27r-periodic and C k functions v : R — > R. Note that the map N maps X into itself. 
We consider the closed subspace of X 

x h = {vex, C» = o}, 

where 

TT J — 7r 

Note that C*(u) = f j^u^) sin^d^ for all u G X. We have X = R sin£ © X fe and 
note Pu = C*i u ) sin ^, = I — P the corresponding projectors. 

Now we split small amplitude solutions as in section [2] using the ansatz ( [16]) and 
the simplified expression (I3ip. in which we fix the phase <fi = —n/2 to recoved odd 
solutions. More precisely, we set 

u(0 = asm£ + a a v(0, (45) 

where a > is assumed close to and v G D h = X h fl Cp er (0, 2ir). This splitting is 
reminiscent of the classical Lyapunov-Schmidt technique, since the closed operator 
at the left side of (14*3]) has a kernel spanned by sin £ in X and its restriction to X^ is 
invertible. The nonlinear dispersion relation (1321) suggests the scaling u 2 — 1 = a a_1 A, 
where A G R and v G Dh will be subsequently determined as functions of a. Inserting 
(1451) in (1431) and projecting the resulting equation on the first Fourier harmonic 
results in 

\ = F (v,a), (46) 
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we project (jl3"|) onto X h , apply K h — (-^ + 1) 1 G C{X h) D h ) to both sides of 
and use the fact that K h -^L = S^K h = I — K h on D h . This yields 



F (v,a) = -a- a C[N{a sin£ + a Q v)]. (47) 
Now there remains to complete ( )46l) by a fixed point equation for v. For this purpose 

j<r, — 

v = F h (v,X,a), (48) 

F h (v, A, a) = a*' 1 X(K h - I)v + a~ a K h P h N(a sin£ + a a v), (49) 
where Kh is explicitly given by equation (jUJ) (checking that Kh maps Xh into D/j 
is lengthy but straightforward). As a conclusion, nontrivial solutions of ( 142|) taking 
the form ( )45i) are given by the fixed point equation 

(v,X)=F a (v,X), (50) 

where the map F a = (F h , F ) is defined by fj47l) - fj49l and depends on the parameter 
a > 0. We shall consider fl50|) as a fixed point equation in some closed ball of Xh x R. 
Then any solution of fl50|) will have the automatic smoothness v G Dh since 



2 - - K h ( a"" 1 Aw + cT Q P fe iV(a sin ^ + a a v) ) (51) 



a; v 



and K h G C{X h ,D h ). 



3.2. Solution of the fixed point problem. In what follows we prove some key 
properties of F a required to apply the contraction mapping theorem. We shall 
note V = (v,X) G Xh x R and equip Xh X M. with the supremum norm || V || co = 
Max(||i;||L«>, |A|). We denote by B(p) the closed ball of center and radius p in 
X h x R. 

Lemma 1. Consider the map F a = (Fh(. , a) , F (. , a)) defined by expressions ffl 7 ?])- 
fijm , where the nonlinear term N takes the form fi44\ ). Assume that the potential V 
defining N satisfies the same conditions as in theorem^ Then there exist p > 0, 
a > independent of q such that F a maps B(p) into itself for all a G (0, do). 
Moreover, one has as a — » + 

sup IIW)-W)IU =OK -. ) 

v 1 ,v 2 eB( P ),v 1 ^V2 ||Vi — mioo 
uniformly in q G [— 7r,7r]. 

Proof. We first prove that _F a maps a suitable ball -B(p) into itself. For notational 
simplicity we shall use the same symbol C for different multiplicative constants that 
are all independent of p, q and a. The assumptions made on V in theorem [1] imply 
property (fT7|) . from which it follows that = O(||u||"oo) as — > 0, 

and thus ||iV(u)||i/x> = 0(||m||£oo) uniformly in q G [— 7r,7r]. Consequently, given 

2 

p > 1 and ao(p) = p 1 ^, there exists C > such that for all a G (0, ao(p)), for all 
(v, A) G B(p), for all g G [— n, it), 

\\a- a N(a sin ^ + a a v)\\ L o C <C, (52) 

\\a a ~ l X(K h - I)v\\ L ~ < C. (53) 
This yields immediately the same type of estimate for F a , namely 

\\F a {v, A) ||oo < C. (54) 
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Since C is independent of p, F a maps B(p) into itself provided p > Max(C, 1) and 
a < a (p). 

Now we fix a value of p such that F a maps B(p) into itself for small enough a and 
estimate the Lipschitz constant of F a . In what follows we use the same notation 
C for different multiplicative constants independent of q and a. The assumptions 
made on V in theorem [TJ and the mean value inequality imply 

\V'{x)-V'{y)\<C\x-y\\\{x,y)\\Z l 

for all (x, y) G M 2 close to 0. Consequently, one obtains successively for all functions 
U!,u 2 G L°°(R) close to 

WV'M - V'(u 2 )\\ L °c < C |K - n 2 || LOO ||(« ls « 2 )||^ )a , 

\\N(ui) - N(u 2 )\\ L oo < C \\u x - n 2 || Loo \\(u u u 2 )\\«~y )2 

uniformly in q G [— 7r, tt]. Consequently, there exists C > such that for all a small 
enough, for all {v j, Aj) G -B(p), 



a 



(iV(a sin^ + a^i) - iV(a sin^ + a a v 2 ) ) < Cfl"" 1 |K - v 2 ||l», (55) 

||F a ( Vlj AO -F> 2 ,A 2 )|U < CTa"- 1 ||(«i, Ai) - (v 2 ,A 2 )|U. (56) 

□ 

As a corollary of lemma HJ one obtains the following result. 

Theorem 2. Fix B(p) as in lemmaUl There exists ao > independent of q such 
that F a admits a unique fixed point V(a) = (v a ,X a ) in B(p) for all a G (0, ao). 
Moreover, there exist M,C > independent of q such that the sequence {V^(a)} n >o 
in B(p) defined by V n (a) = F" +1 (0) satisfies 

\\V n (a) - Via)^ < MiCa"- 1 )^ 1 for all n > 0. (57) 

Proof. According to the estimate of lemma HJ F a defines a contraction on B(p) for 
all a small enough, with Lipschitz constant k < C a Q_1 , where C is independent of 
q. In that case, by the contraction mapping theorem, F a admits a unique fixed point 
V(a) in B(p) which is exponentially attracting. Moreover, the classical estimate of 
the convergence speed yields 

n+l 

||K(a)-^(a)|U<:j -IIWIU 

which implies ( 157|) for a m 0. 

□ 

Theorem [2] implies the existence of nontrivial solutions (u a ,ul) of (1421) taking the 
form 

u a (0 =a sm£ + a a v a (0, (58) 

u 2 a = 1 + a a ~ l X a . (59) 

This in turn implies the existence of periodic travelling wave solutions of ([T]) taking 
the form (j3j) . In order to complete the proof of theorem [TJ we now need to determine 
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the principal parts of v a and A a for a ~ 0. Considering F a (0) = (t> (a), Ao(a)), we 
have by estimate fl57|) 

\\v {a) - v a \\ L o. + |A (a) - A a | = 0(a a - r ), (60) 

hence we are led to compute vo(a) and Ao(a). One finds 

Ao(a) = -a- Q C*[AT(asinO], (61) 

v (a) = a- a K h P h N(a sin f), (62) 

where Kh is defined by ([9]). Keeping in mind expansion (ITT)) of V, the principal 
part of the nonlinear map N reads 

[No(u)](0 = VMt + - VMO - «(£ - <?)), (63) 

where V^'(x) = — -H"(— x). For all a > one has V^'(ax) = a a V^(x), and in the 
same way No(au) = a a No(u). Using this property and expansion (ITT)) in equations 
()6Tl) - (l62l) . one obtains consequently as a — > 

A (a) = -C[iVo(sinO]+o(l), (64) 
v (a) = K h P h N (smO + o(l) in £> h . (65) 

Estimate (160)) and expansion ( 165)) yield expansion (JTJ of theorem [TJ Now there 
remains to recover expansion (jSJ), which requires to compute leading order terms 
of (l64p . We can restrict our computations to the case q G [0, 7r], which allows one 
to recover the results for q G [— 7T, 0] by symmetry considerations (see section 133)) . 
After lengthy but straightforward computations based on classical trigonometric 
identities, one obtains for all q G [0, n] 

C[^(sin(e + g)-sinO] 
= -2[2sin(^r f |cos(£ + ^)| Q sin^ 

7T Z I TV q 2 

= -- [2sin (§)]"( sin (§) f +2 | cos^ 1 ds + -^[cos (fr+ 2 ), 



7T V 2 /J v v 2 7 J* 1 a + 1 L v 2 

2 

and in the same way 

CTC(sin£-sm(£-g))] 

7T 

= -[2sin(fr(sin(|) / 2 (cos,r +1 rf, - ^-[cos (^)] Q+2 ). 

7T 2 2 1 3. OL + 1 2 

2 

Combining (l63l)-( l64l) with the above identities, one finds after elementary computa- 
tions 

A (a) = [2sin(|)] a+1 c a + o(l), (66) 

where c a is defined in equation (J2TJ). Combining expression ( 159)) with estimate (160)) 
and expansion (166)) yields finally for q G [0, 7r] 

w 2 = l + a Q - 1 [2sin(|)] cl+1 c Q + o(a cl ~ 1 ). (67) 



Expressing c a using Euler's Gamma function as in (f23j) . one finally obtains expansion 
(JSJ) of theorem UJ for g G [0, 7r]. 
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Remark 4. As shown by estimate (57]), computing V n (a) forn large enough allows 
one to obtain approximations of the fixed point (v a , X a ) at any order in a m 0, which 
provides approximations of the wave profile and wave frequency at arbitrary or- 
der. However the expression ofV n (a) becomes highly complex already at n = 2, hence 
this series of explicit approximations doesn't seem useful in practice. Nevertheless, 
it is interesting to see that approximations at arbitrary order in a can be obtained 
despite the limited smoothness of V at the origin. 

Remark 5. It is again interesting at this stage to reinterpret the results in the 
framework of Lyapunov-Schmidt reduction. Adding an arbitrary phase <fi to solu- 
tion and introducing A = ae 1 ^, the nonlinear dispersion relation (5)) describes 
nontrivial solutions of the bifurcation equation 

A (u - 1) A {Al*- 1 1 sin (-)| a+1 + h.o.t. = 0. (68) 

T 2 

More precisely, solutions (fJ])-([3J) and their phase shifts appear through a pitchfork 
bifurcation with 5*0(2) symmetry which occurs at (A,oj) = (0,1). Note that the 
classical pitchfork bifurcation equation corresponds to the case a = 3 of Iffify) , and 
for a > 1 it can be recovered from ( ESP after the C° change of variable Z = \A\^~ A. 

3.3. Symmetries. In this section we point out some symmetry properties of the 
periodic travelling waves. We begin by examining more closely their dependency 
on the wavenumber q, hence we will denote the map (|44p by N q , the map F a by 
F a>q and the fixed point of F a>q in B(p) by V(a,q) = (v atq ,X atq ). Let us introduce 
the operator [r vr M](^) = u(£ + tt) corresponding to half-period space-shift and the 
symmetry S G C{X) defined by S = —r n . One can check that 

SN q (a sin^ + a a v) = N_ q (a sin£ + a a Sv). (69) 

Consequently, from the identity (v a>q , A 0i9 ) = F a q {y a ^,\ a ^ q ) and definition (j4"B1) - (l4"g]) 
of F a we obtain 

(S v a>q , X a>q ) = F a _ q (S v a>q , X a:q ) 

since S commutes with and P S = S P = P. By the uniqueness of the fixed 
point of F a - q in B(p), it follows 

v a - q = S v ayq , X a _ 9 = X a>q . (70) 

In particular, expression ( 1671) becomes for all q G [— tt, tt] 

u 2 a = l + a a ^\2 sin (|)| Q+1 c a + o(a a - v ), (71) 

hence one finally recovers expansion (jSJ) of theorem [1] for all q G [— 7T, tt]. 

In what follows we concentrate on the particular case q = it. Since functions 
belonging to X are 27r-periodic, the operators N q and F a>q (acting respectively on X 
and X/jXR) are 27r-periodic in q, and the same holds true for v a>q , X a>q . Consequently, 
property (1701) yields v a>n = Sv a>w , i.e. for all £ G R one has 



(72) 
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It follows that for q = it, the solution (jl]) takes the form 

x n (t) = a (-l) n+1 (sin (cu a t) + a^v^uj)), (73) 

i.e. it corresponds to a binary oscillation for which nearest-neighbours oscillate out 
of phase with opposite displacements. 

Note that there exists another simple method to obtain binary oscillations in the 
special case V — Vq, without restriction to the small amplitude limit. One can find 
solutions of ([TJ) of the form x n (t) = (— l) n a(t) after a short discussion with respect 
to the sign of a(t) and parity of n. This expression defines a solution of ([T]) if and 
only if a satisfies the integrable equation 

^ + a + 2 a a\a\ a - 1 = 0, (74) 
at 2 

the phase space of which is filled by periodic orbits. 

4. Numerical comparison of the DpS model and Newton's cradle 

In this section we numerically analyze how the dynamics of the DpS equation 
approximates the one of the nonlinear chain ([1]), in the case of classical Hertz's 
contact law (V — V , a — 3/2) and for some classes of initial conditions. We 
start our investigations with initial conditions leading to periodic travelling waves. 
In section [3] we have proved the existence of exact periodic travelling waves of ([T]) 
whose profile and velocity are determined at leading order by the DpS equation in 
the small amplitude limit. However these results do not indicate the range of validity 
of the asymptotic analysis and do not provide informations on the local dynamics 
around the travelling waves. We examine such questions in section 14.11 for stable 
travelling waves and in section [4T21 for unstable ones. As we shall see, system ([1]) and 
the DpS equation display similar features concerning the occurence of modulational 
instability, which yields the formation of discrete breathers. A more precise study of 
these localized structures requires well chosen localized initial conditions. In section 
14.3} we use initial conditions deduced from the study of localized standing waves of 
the DpS equation (section 12.21) in order to generate static and travelling breathers. 

Let us now describe our general numerical procedure before going further. We 
consider a given initial condition {Ai(0)} n gz and use ansatz ffTTj) to determine an 
initial condition for ([[]) (computing x n (0) requires the knowledge of d T A n (0) which is 
determined by the case r = of ffTUl) ). Then we integrate ([T]) and (fTUj) numerically 
(or use an explicit solution of ffTU]) when available) and compare the solution of 
(JTJ) with the ansatz (|TT|) at subsequent times. In what follows the corresponding 
solution of (|TJ) will be referred to as the exact solution, and the approximate solution 
will denote the ansatz (|TT|) computed with the DpS equation. 

We realize different tests for chains ranging from 50 or 200 beads with periodic 
boundary conditions. Numerical integrations are performed using the standard ode 
solver of the software package Scilab. The numerical integrator doesn't exactly 
conserve the Hamiltonian ([3]), but in all simulations we have checked that the relative 
errors \H(t) — 7^(0)1/^(0) are less than 3,5. 10~ 4 , this upper bound corresponding 
to simulations over long times of the order of 540 periods of the linear oscillators 
(section 14.3ft . 
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4.1. Approximation of stable travelling waves. We begin by computing peri- 
odic travelling waves in the small amplitude regime in which the weakly nonlinear 
analysis of section 12.11 has been performed. We normalize the travelling wave solu- 
tion ( 128]) by fixing R — 1, so that a = 2e in ansatz (13T1) . For the first numerical run, 
we choose (131]) as an initial condition, with <f> = 0, q = n/5 and u tw = 1.01, which 
yields a ~ 2.10 -2 in the nonlinear dispersion relation (I3"2"j) . These parameter values 
correspond to a weakly- nonlinear regime, since the maximal value of the nonlinear 
interaction forces (V^x^ — x^)\ estimated with (T5T|) is close to a 3 / 2 /2, which cor- 
responds to seven percent of the maximal linear restoring force a. As shown by 
figure [5] (left plot), the agreement between the exact and approximate solutions is 
excellent, even over long times corresponding to 100 multiples of the wave period 
T tw = 2n/uJt w . More precisely, the relative error between exact and approximate 
solutions (measured with the supremum norm) is less than five percent at the final 
time of computation (figure 0, right plot). These results show that the wave profile 
remains almost sinusoidal. However, nonlinear effects cannot be neglected on the 
timescales we have considered, because they still affect the wave velocities. This ap- 
pears clearly in the right plot of figure EJ which compares the exact and approximate 
solutions to the linear wave with same amplitude a and frequency to — 1 < u tw . This 
wave has a lower velocity than the exact solution, so that it becomes approximately 
out of phase at time to = 333.7. Nonlinear selection of the wave velocity is cor- 
rectly captured by the DpS equation, since the graphs of the exact and approximate 
solutions are almost superposed. 

In a second numerical run, we modify the previous initial condition by setting 

A„(0) = (1 + pW) cos (qn) - i (1 + p«) sin (qn), (75) 

where pn,Pn £ [— 1/2, 1/2] are uniformly distributed random variables. The mis- 
match between exact and approximate solutions increases significantly compared to 
figure [5] (see figure El left plot), but the approximation remains very satisfactory 
given the large time interval considered. A better agreement between the exact and 
approximate solutions could probably be achieved by taking into account the higher 
order correction R n to ffTTT) provided by f|T8|) . Besides the comparison between exact 
and approximate solutions, this numerical experiment reveals the stability of the 
corresponding travelling waves of equations (pQ) and ffTUl) . at least over 100 periods 
of the linear oscillators (however instabilities might occur on longer timescales). As 
we shall see later, fast modulational instabilities can show up for other choices of 
wavenumbers q. 

Next we again choose ansatz (I3T1) with q = it/ '5 as an initial condition, but 
increase the frequency up to u tw = 1.1, which yields a ~ 2.11 in (I321) . In that 
case, the maximal values of the nonlinear interaction force and the linear restoring 
force estimated with (I3TI) are of the same order, i.e. this initial condition should 
correspond to a fully nonlinear regime. As shown by figure [7J the initial condition 
still generates a periodic travelling wave that is qualitatively well described by the 
DpS equation. Quite surprisingly, the waveform remains nearly sinusoidal and close 
to the ansatz (13 ip (see figure [S]). However the ansatz doesn't accurately describe the 
wave velocity, as revealed by the shift between the exact and approximate solutions. 
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Figure 5. Left plot : small amplitude exact solution (continuous 
line) and approximate solution (dash-dot line) for an initial condition 
corresponding to ansatz (pTTj) . with = O,g = 7r/5 and u tw = 1.01. 
Bead displacements are plotted at a time t = 621.95 of the order of 
100 periods of the linear oscillators. The agreement is excellent since 
both graphs are almost indistinguishable. Right plot : relative error 
E(t) = a" 1 \\x^(t) — (t) 1 1 oo between exact and approximate solu- 
tions. 




Figure 6. In the left plot, the same experiment as in figure |5] is per- 
formed with a random noise added to the initial condition (equation 
(1751) ). The right plot compares the exact and approximate solutions 
to the linear wave x n (t) = a cos (qn — t) solution of equation (JTJ lin- 
earized at x n = 0. Comparison is made at time to = 333.7. This wave 
has a nearly identical profile (dashed line) but a lower velocity, so 
that it becomes out of phase with the exact solution at t — to- O n the 
contrary, the approximate solution deduced from the nonlinear DpS 
equation captures the correct wave velocity with excellent precision, 
and the graphs of the exact and approximate solutions (continuous 
and dash-dot lines) are almost superposed. 

In the next computation we keep the same parameter values except we increase the 
frequency up to u tw = 1.5, so that a 52.9. This correspond to a highly nonlinear 
regime, where one can expect that the DpS equation will not properly describe the 
dynamics of (pQ). The two plots of figure [9] reveal indeed notable differences between 
the travelling wave patterns of the exact and approximate solutions. The initial 
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Figure 7. Space-time diagram showing the interaction forces 
Vo(x n +i — x n ), for an initial condition given by ansatz (1ST]) with 
q = ii / 5 and co tw = 1.1. Forces are represented in grey levels, white 
corresponding to vanishing interactions (i.e. beads not in contact) 
and black to a minimal negative value of the contact force. The left 
plot corresponds to the exact solution, taking the form of a periodic 
travelling wave with almost constant velocity. The right plot shows 
the approximate solution. 




Figure 8. Wave profiles of the exact (continuous line) and approx- 
imate (dash-dot line) solutions in the numerical experiment of figure 
0, shown at two different times t = 64.95 (left) and t = 570.95 (right). 
For the exact solution, the travelling wave motion is coupled to small 
collective in-phase oscillations, which explains the vertical shift of the 
profile. The wave velocity is less accurately described by the DpS 
equation in this regime, so that a small shift of the two profiles is al- 
ready present at t = 64.95, and the exact and approximate solutions 
are out of phase at t — 570.95. 



condition generates in fact a pulsating travelling wave, showing internal oscillations 
that appear as spots along the wavefronts of figures 191 and [TUl This internal breathing 
appears at the beginning of the simulation (around t = 1) and releases some energy 
in the form of small amplitude travelling waves with lower velocity. These waves 
subsequently collide with the main travelling wave, generating travelling waves with 
negative velocities. These phenomena are absent in the DpS equation as shown by 
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the right plot of figure [9j In addition the wave profile is not any more sinusoidal 
(see figure [TUl right plot). 




FIGURE 9. Same computation as in figure [3, but this time in a highly 
nonlinear regime where the frequency parameter of the initial ansatz 
is increased to u tw = 1.5. In that case the travelling wave pattern of 
the exact solution (left plot) develops a complex small-scale structure 
that is absent in the approximate solution (right plot). 
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Figure 10. Left plot : zoom on the left-side space-time diagram of 
figure showing the small-scale structure superposed on the travelling 
wave pattern of the exact solution. Right plot : comparison between 
bead displacements for the exact (continuous line) and approximate 
solution (dots) at t = 83.7. 

As a conclusion, we have seen that the initial condition obtained with ansatz 
(13"TT) and q = n/5 generates a stable travelling wave that is qualitatively well de- 
scribed by the DpS equation if one avoids the highly nonlinear regime. Moreover, 
in the weakly nonlinear regime the DpS equation describes the travelling wave with 
excellent precision. 

4.2. Modulational instability and breathers. In what follows we perform the 
same kind of simulations but increase the wavenumber q above tt/2, which yields 
a completely different dynamical behaviour corresponding to modulational instabil- 
ity. We first determine an unperturbed initial condition using ansatz (I31I) - (I32I) with 
q = Att/5 and u tw = 1.1, which corresponds to a small amplitude e ~ 3,8. 10~ 3 in 
(1301 (as above the travelling wave solution ff28l) is normalized by fixing R — 1). Then 
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we perturb this initial condition as indicated in ( 175]) . where pn\pn £ [—0.1,0.1] 
denote uniformly distributed random variables. The phenomenon of modulational 
instability is illustrated by figures [11] and [12] The envelope of the initial condition 
(figure [11] left plot) localizes after some transcient time, yielding larger oscillations 
of some beads corresponding to discrete breathers (figure [TT] right plot). The mod- 
ulational instability yields a disordered train of travelling breathers that propagate 
along the lattice as shown by figure [TU The DpS equation describes this modula- 
tional instability quite well (compare the two plots of figure [TJ]), and even reproduces 
fine details like the crossing of two breathers at t « 330. 




Figure 11. Left plot : initial bead displacements, fixed by the pe- 
riodic travelling wave ansatz (|3~T]) -( 13"2"]) with q = 47r/5 and u tw = 1.1, 
perturbed by a small random noise. Right plot : bead displacements 
at t = 331.2. The dynamics generates alternate regions of large and 
small amplitude motion, characteristic of modulational instability. 



Figure 12. Left plot : spatiotemporal evolution of the interaction 
forces for the computation of figure [TTJ Forces are represented in grey 
levels as in figure [7J The grey region at early times of the simulation 
corresponds to a travelling wave profile, and the dark lines to a train of 
travelling breathers generated by the modulational instability. Right 
plot : approximate solution deduced from the DpS equation. 



The average velocity of the travelling breathers generated near the onset of insta- 
bility depends on the wavenumber of the unstable travelling wave. When q is equal 
or close to n (q = n corresponding to binary oscillations), one typically observes 
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static breathers that remain pinned to some lattice sites over long times. This phe- 
nomenon is illustrated by figures [13] and [TU which correspond to the case q = ir, 
uj tw = 1.1 and ph 1 iPn £ [—0.5,0.5]. The DpS equation provides a good qualitative 
picture of the instability (compare the two plots of figure U4|) . It is able to reproduce 
some complex phenomena like the interaction and subsequent merging of two close 
static breathers, at t ~ 223. However the local dynamics after the two breathers 
have merged is not well described by the DpS equation, since more energy is released 
at the collision site for the exact solution than with the approximate solution. This 
may be the sign that higher harmonics neglected by ansatz (11 II) come into play dur- 
ing the merging phase, or may simply originate from sensitiveness of the collision 
result to relative phases and breather positions. 

In the above space-time diagrams, travelling and static breathers appear as grey 
lines due to their internal breathing motion. These oscillations are more visible on 
figure [TBI where static breathers with different amplitudes and spatial extentions 
are shown. 




n n 



Figure 13. Left plot : initial bead displacements, corresponding to a 
binary oscillation (q = ir and u tw = 1.1 in ansatz fl3~Tj) - fl32"]) ). perturbed 
by a random noise. Right plot : snapshot of a static breather formed 
at t ~ 223 after the merging of two smaller breathers. 



During the modulational instability, only a part of the modes are initially ampli- 
fied. In order to analyze how the DpS equation is able to reproduce this phenom- 
enon, we now study the time evolution of the spatial Fourier transform of {x n (t)}. 
As previously we consider periodic boundary conditions x n+ ^(t) = x n (t), but in- 
crease the number of particles to iV = 200 in order to achieve a better spectral 
resolution. Considering the discretization of the spectral band [—ir, vr] given by 
= • { — y + l,...,y} and defining the discrete Fourier transform 

N-l 

x k (t) = Y,Xn(t)e- ink , fcer„ 

n=0 

we have 
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Figure 14. Spatiotemporal evolution of the interaction forces for the 
computation of figure [T51 (left plot), and its approximation deduced 
from the DpS equation (right plot). The grey region at early times of 
the simulation corresponds to the perturbed binary oscillation, and the 
dark regions to nearly static breathers generated by the modulational 
instability. 



Contact forces 



Figure 15. Detail of figure HH showing a distribution of static 
breathers with different sizes and amplitudes. Black spots alternat- 
ing with white regions correspond to sequences of compressions and 
separations of neighbouring beads. 

We compare the time evolution of the magnitude of the Fourier transform \xk(t) | for 
the exact and approximate solutions, for initial travelling wave profiles perturbed 
by random noise. The unperturbed initial condition is determined using ansatz 
fl3ll - fl32l) . for (jj tw = 1.1 and different wavenumbers q that determine the amplitude 
e = a/2 in fl30|) . We perturb this initial condition as indicated in f J75|) . where 
Pn j pff are random variables generated from the Gaussian distribution with mean 
and standard deviation 0,5. 10~ 2 . 

Figure [16] compares the spectra of the exact and approximate solutions for q = 
47r/5. At t = 106.7 (first row), the two peaks at k = ±q correspond to the periodic 
travelling wave generated by the initial condition. Two second harmonics at k = 
±2(g — 7r) are visible for the exact solution and are absent for the approximate 
one (these harmonics are generated at early times of the simulation even in the 
absence of noise). As time further increases, a band of unstable modes grows in 
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amplitude at both sides of k = ±q. The unstable bandwidths and most unstable 
wavenumbers of the exact and approximate solutions are quite close (compare the 
plots at t — 399.95, third row), but the approximate solution overestimates the 
initial growth rate (compare the plots at t — 253.45, second row). Both for the 
exact and approximate solutions, a breathing of the spectrum near k = ±7r settles 
progressively. The plots at the fourth row show quite similar spectral distributions 
of the exact solution at t = 541.95 and of the approximate one at t = 543.7 (we have 
slightly shifted time to compensate a phase-shift in the spectral breathing). As a 
conclusion, the DpS equation describes quite accurately the modulational instability 
of the travelling wave with q = 4tt/5 in the weakly nonlinear regime. 

In contrast with the above results, we describe below two situations in which 
the DpS equation gives incorrect results on the stability or dynamical properties of 
travelling waves. 

The first one corresponds to a simple instability threshold effect. For q = 7r/2 and 
below, we haven't observed modulational instabilities both for the DpS equation 
and system ([T]) in the weakly nonlinear regime. However, in both models we have 
found thresholds for the modulational instability at slightly different critical values 
of q near ir/2. This is illustrated by figure [13 where the magnitudes of the spatial 
Fourier transforms of the exact and approximate solutions are plotted at t = 913.7 
(160 time periods of the initial periodic travelling wave) for different values of q near 
7r/2. As shown by the two plots at the first row, small amplitude travelling waves 
with q = 527r/100 are unstable both for system ([1]) and the DpS approximation. 
When q decreases to q = 5l7r/100, the modulational instability persists for the DpS 
approximation, but it does not occur (or may occur extremely slowly) for system 
([TJ) (see plots at the second row). The third row of figure [T71 provides the spectra 
for q = 7r/2, showing no trace of modulational instability. In addition to the main 
peaks, one can see second harmonics at k = ±ir in the exact solution spectrum, and 
a small noisy background (more important for the exact solution) that follows from 
an energy cascade between modes. As a conclusion, these computations reveal that 
system (Cp) and the DpS approximation yield slightly different instability thresholds 
for the wavenumber q, and consequently the DpS approximation does not correctly 
describe the dynamics when q is chosen between the two critical values. 

As second situation in which the DpS approximation breaks down is described in 
figure [18j For q = 577r/100, the travelling wave develops a modulational instability 
both in system ([1]) and the DpS equation, but the two models display different tran- 
sitory dynamical behaviours. This originates from two small harmonics appearing in 
the spectrum of the exact solution at early times of the simulation, more precisely a 
second harmonic at k = ±2(g — tt) and a third harmonic at k = ±(3g — 2n). As pre- 
viously these peaks are absent in the spectrum of the approximate solution. Here the 
difference with respect to figures [H] and [TT] is that these small peaks fall inside the 
bands of modes amplified by modulational instability, and their amplitude grows 
with time (the second harmonic dominates the third one at the beginning of the 
simulation, but their amplitudes become comparable arount t = 150). This results 
in the spectra shown in figure [TH where two additional large amplitude oscillating 
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Figure 16. Magnitude of the discrete spatial Fourier transforms 
of the exact solution (left column) and approximate solution (right 
column), for a random perturbation of a small amplitude travelling 
wave with q = 4-7r/5. Plots correspond to different times t = 106.7 
(first row), t = 253.45 (second row), t = 399.95 (third row), t = 541.95 
(fourth row, left), t = 543.7 (fourth row, right). 



s are present in the spectrum of the exact solution at t — 344.95 (upper left 
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Figure 17. Magnitude of the discrete spatial Fourier transforms of 
the exact solution (left column) and approximate solution (right col- 
umn), for a random perturbation of small amplitude travelling waves 
with different wavenumbers q close to 7r/2 (top plots : q = 527r/100, 
middle : q = 5l7r/100, bottom q = ir/2). Plots correspond to time 
t = 913.7. 



and absent for the approximate solution (upper right plot). The two middle plots 
show the bead displacements for the exact solution at two different times, revealing 
a large scale standing wave that strongly modulates the basic pattern. The lower 
plot corresponds to the approximate solution given by the DpS equation, which is 
very weakly modulated and fails to reproduce the correct transitory dynamics of ([T]) 
in the present situation. However, in both models broad bands of unstable modes 
become slowly amplified at larger times by modulational instability (this appears 
clearly after t = 600 for system (pQ) and t = 700 for the DpS equation), and around 
t = 900 the bead displacements show rather similar features in both cases. 
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Figure 18. Upper plots : magnitude of the discrete spatial Fourier 
transforms of the exact solution (left column) and approximate solu- 
tion (right column), for a random perturbation of a small amplitude 
travelling wave with wavenumber q = 577r/100, at time t = 344.95. 
Middle plots : bead displacements for the exact solution at t = 402.95 
and t = 407.2. Lower plot : bead displacements for the approximate 
solution at t = 407.2. 



4.3. Localized initial conditions. In section 14721 we have seen that the modula- 
tional instability of certain travelling waves generates travelling or static breathers, 
a phenomenon which can be captured by the DpS equation. In what follows we 
illustrate with a few examples how the DpS equation approximates the evolution 
and interaction of these localized structures. 

We start by generating single breathers similar to the ones that result from the 
modulational instability of the q = n mode. For this purpose, we consider the 
spatially antisymmetric homoclinic solution {a n } of fl37|) represented in figure H] 
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(bottom left plot), the associated breather solution of ([24"]) 

Mr) = \^a n e l \ fi = (u sw - 1) 2% > 0, (76) 

and the corresponding breather ansatz (HIT) . Using the initial condition determined 
by (jlTjl . we integrate ([1]) numerically. Computations are performed for a chain of 
iV = 50 particles with periodic boundary conditions. 

The results are shown in figure [19] for different amplitudes of the initial breather 
profile obtained by varying /i. The first case corresponds to |xjv/2+i(0)| ~ 9,7.10~ 3 , 
i.e. u sw ~ 1.1 (initial condition shown in the top left plot). By construction, the 
profile of the approximate solution given by the DpS equation remains similar to 
the top left plot for all times, just oscillating periodically in time. One observes the 
same qualitative behaviour for the exact solution. As shown by the top right plot, 
the profiles of the exact solution at t = 569.4 (continuous line) and the approximate 
solution at t = 565.95 (dots) are very close (we compare both profiles with a small 
time-shift, because the exact and approximate solutions develop a half-period phase- 
shift after approximately 83 breather periods). These results show the efficiency 
of approximation f ]4"T]) on large time intervals of the order of 100 multiples of the 
breather period. As shown by the bottom left plot, this localized 3-sites breather 
solution qualitatively corresponds to the central breather forming in figure \T5\ In 
the second case, the initial profile has a larger amplitude |xjv/2+i(0)| ~ 0.34, which 
corresponds to u sw ~ 1.6. The bottom right plot (bead displacements at t = 293.2) 
show that small dispersive waves are emitted by the exact solution. On the contrary, 
the DpS equation does not yield any dispersion in that case, only time-periodic 
oscillations of the profile given by the initial condition. 

Remark 6. We obtain qualitatively similar results with the spatially symmetric 
homoclinic solution {a n } of (37\ ) (figure^ top left plot), except the exact solution 
rapidly develops a slight asymmetry which is absent in the approximate solution. 
This is due to the fact that (QP does not possess the invariance x n — > X- n , whereas 
nD\) has the invariance A n — > A_ n . The asymmetry of the approximate solution 
could be recovered by taking into account the higher order term R n present in [W\l . 

The above computation shows that the DpS equation can approximate single 
breather solutions of ([T]) provided their amplitude is sufficiently small, i.e. their 
frequency sufficiently close to 1. Now let us illustrate how the DpS equation accounts 
for the interaction of small amplitude breathers, on long times scales up to 600 
breather periods. Initial conditions for ([1]) are generated using the ansatz ( 14T|) with 
u sw ~ 1.1, for initial envelopes {a n } that do not correspond to exact solutions of 
equation ( 137]) . The DpS equation is initialized by setting r = in ( ITS"]) . 

The left side plots of figure [20] correspond to initial displacements almost com- 
pactly supported on 12 lattice sites, modulating binary oscillation with amplitude 
close to 0.016 (top plot). Initial velocities are set to 0. The spatiotemporal evolu- 
tion of the interaction forces is shown for equation ([T]) (middle plot) and the DpS 
approximation deduced from ansatz ( II ip (bottom plot). In both cases, the initial 
packet splits into a static breather and two pairs of breathers travelling slowly in 
opposite directions. Further crossings of the pulses occur, with some interactions 
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Figure 19. Top left plot : initial condition x n (0) determined by (|4lT) . 
where the amplitudes a n correspond to a spatially localized solution 
of fl37j) and initial velocities are set to 0. This initial condition defines 
an approximate breather solution with frequency u sw ~ 1.1. Top 
right plot : comparison of the exact solution at t — 569.4 (continuous 
line) with the approximate solution at t — 565.95 (dots). Bottom left 
plot : spatiotemporal evolution of the interaction forces for the exact 
solution (grey levels), where a succession of black spots reveals the 
breathing dynamics. Bottom right : exact solution at t = 293.2, for 
the same type of initial condition with amplitude approximately 30 
times larger than in the previous case (u sw ~ 1.6). 

resulting in small phase shifts or velocity changes. The global sequence of collisions 
is different for the exact and approximate solutions (because collisions result in dif- 
ferences that accumulate with time), but these computations show nevertheless a 
good qualitative agreement between the exact and approximate dynamics. 

The right side plots of figure [20] are obtained for initial displacements correspond- 
ing to two small amplitude static breathers initially separated by one lattice site 
(top plot), with amplitude close to 0.01 and zero initial velocities. Merging of the 
two breathers occurs after a long transcient, both for the exact (middle plot) and 
approximate (bottom plot) solutions. However this occurs at largely different times, 
around 175 multiples of the breather period for the exact solution, and twice more 
for the approximate solution. In that case, the DpS approximation gives conse- 
quently a good qualitative picture of the slow interaction of the two breathers, but 
is quantitatively not satisfactory to predict the finite lifetime of the two-breather 
bound state. 
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Figure 20. First row : bead displacements corresponding to two 
different initial conditions (initial velocities are set to 0). Second row : 
spatiotemporal evolution of the interaction forces of system (JTJ for the 
above initial conditions. Third row : interaction forces corresponding 
to the approximate solutions deduced from the DpS equation for each 
initial condition. 

Let us end this section with some remarks concerning boundary conditions. So far 
we have computed discrete breathers in chains of N beads with N > 50 and periodic 
boundary conditions. From an experimental point of view, system (TfJ) corresponds 
in that case to a closed ring of N beads, where N must be large enough in order 
to avoid curvature effects. However most experiments with Newton's cradles have 
been performed with a small number of beads arranged linearly, which corresponds 
to system (TTJ) with free end boundary conditions. It is important to stress that 
discrete breathers can be also generated in this context, as shown by figure I2T1 The 
left plots show displacements in a chain of 10 beads at different times, obtained by 
numerical integration of system ([T]) with free end boundary conditions. The results 
show the existence of a spatially antisymmetric breather involving mainly 4 beads. 
The initial condition is generated using the same method and parameter values as 
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for the small amplitude breather of figure [19] (top left plot), except we use free end 
boundary conditions in the stationary DpS equation ([37]) to compute the breather 
ansatz fj4TT) . The right plots show the approximate solution (jUJ) deduced from the 
DpS equation. The latter agrees very well with the exact solution, except both differ 
by a phase-shift that is increasing with time. 





Figure 21. Left plots : bead displacements for the exact solution, 
obtained by integrating ([T|) numerically with 10 beads and free end 
boundary conditions. The initial condition is determined by ( 14"T1) . 
where the amplitudes a n correspond to a spatially localized solution 
of equation (1371) with free end boundary conditions. Right plots : 
approximate breather solution f [4T]) with frequency u sw ~ 1.1 deduced 
from the DpS equation. Beads are labelled with an index n ranging 
from 1 to 10, indicated on the vertical axis. Their displacements are 
plotted in different time intervals, at early times of the simulation (top 
row) and after approximately 250 breather periods (bottom row). The 
scale of bead displacements is indicated at the top of the vertical axis. 



5. Conclusion 

Although Newton's cradle has a long history which dates back to the 17th century 
[40] , it is still considered as a benchmark in recent research at the cutting edge of 
impact mechanics [521 [53]. Moreover, this system turns out to be very interesting 
from the modern perspective of nonlinear science. It allows a simple experimental 
realization of three nontrivial nonlinear effects, namely fully nonlinear dispersion 
[64J , a limited smoothness of nonlinear interactions and their unilaterality. The last 
two aspects are still rather unexplored in the context of nonlinear waves. As we 
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have seen in this paper, these different features require to revisit several classical 
topics, i.e. local bifurcations of periodic waves, modulation equations with fully 
nonlinear dispersion described by discrete p-Laplacians, modulational instability 
and the bifurcation and interaction of discrete breathers. These topics are obviously 
very broad, and our work is only a first step towards a better understanding of these 
questions in a new context. 

From an analytical point of view, we know by theorem [T] the existence of exact 
periodic travelling waves of (DQ) close to small amplitude approximate solutions given 
by the DpS equation. The above numerical simulations have revealed that the DpS 
equation has in fact a much wider applicability. Indeed, it captures other important 
features of the dynamics of (pQ) in the weakly nonlinear regime, namely modula- 
tional instabilities, the existence of static and travelling breathers, and several types 
of interactions of these localized structures. An interesting open question is to jus- 
tify mathematically this close relation between the dynamics of ([T]) and the DpS 
equation, for well-prepared initial conditions close to (ITT]) at t = and large finite 
time intervals. Interestingly, such results have been proved in a different context for 
the DNLS equation (IT21) . where this system is used to approximate solutions of the 
Gross-Pitaevskii equation with a periodic potential (see (62] and references therein). 

In the same spirit, it has been shown recently [33], [3U EJ |66] that the continuum 
cubic nonlinear Schrodinger (NLS) equation (which can be derived in continuum 
limits of the DNLS equation) approximates the evolution of well-prepared initial 
data in many nonlinear lattices. It would be interesting to know if the above results 
extend to equation (|27|) or its generalizations to finite-wavelength continuum limits 
of the DpS equation. Our situation is more complex than in the NLS case, due 
to the limited smoothness of (JTJ) at the origin, and because the Cauchy problem 
for (12 7p seems delicate to analyze. In addition, the analysis of continuum limits of 
the DpS equation may provide interesting informations on the nonlinear waves of 
([1]), e.g. additional explicit approximate solutions. In particular, explicit compactly 
supported travelling or standing waves have been found recently in some variants 
of equation (j27j) that were formally introduced in [7U [75] and in a continuum limit 
of the D-QLS equation [61] (see remark [I] p|8]). Such solutions may also exist for 
continuum limits of the DpS equation, and may accurately approximate the almost 
compactly supported breathers computed in section 12.21 and 14.31 

Interesting open problems related to discrete breathers in Newton's cradle ([I]) and 
the DpS equation concern the analytical proof of their existence, their numerical 
continuation, as well as their stability, movability and interaction properties. It 
will be interesting to compare their collision properties in both models, a problem 
which requires a statistical treatment due to its sensitiveness to relative phases 
and breather positions [9J. In the context of impact mechanics, discrete breathers 
may be used for stringent numerical tests of multiple impact laws in model ([T]). 
They may be also useful for practical purposes, since granular chains are interesting 
devices for the design of shock absorbers (see [611 ES] and references therein), and 
the possibility might exist to trap a part of the energy of an incident wave into static 
breathers in an efficient way. However, before comparing our theoretical results to 
experiments it will be necessary to evaluate the impact of dissipation [12] and spatial 
inhomogeneities on the nonlinear waves we have considered. 
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Another interesting question is the study of the modulational instability of trav- 
elling waves in the DpS equation and its more exhaustive comparison with modu- 
lational instabilities in system ([1]). One of the great interests of the DpS equation 
stems from the fact that the first part of the problem can be worked out analytically, 
as done in reference [15] for the DNLS equation. In addition a more precise numer- 
ical study of modulational instabilities in Newton's cradle and the DpS equation 
(comparing e.g. the most unstable modes in both models) will require statistics on 
power spectra of perturbations. 

More generally, one can wonder if the DpS equation or higher-dimensional ex- 
tensions could capture some general features of nonlinear waves in more general 
granular media. Interesting questions related to nonlinear waves arise in this con- 
text, e.g. the possibility of earthquake triggering by travelling or standing waves in 
granular fault gouges jU]. 
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